package flexibilityDesign;

import ilog.concert.IloConversion;
import ilog.concert.IloLinearNumExpr;
import ilog.concert.IloNumVar;
import ilog.concert.IloNumVarType;
import ilog.concert.IloRange;
import ilog.cplex.IloCplex;
import ilog.cplex.IloCplex.DoubleParam;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.io.Writer;
import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Random;

import concurrent.Executer;
import util.Lang;
import util.NullOutputStream;

public class TestMulticut {
	private int n;// plant number
	private int m;// product number
	private int simulationCount;// sampling number

	private double[][] demands;// clustered demands
	private double[] capacities;// capacities of plants
	private double[][] link_c;// link cost
	private double[][] profit_m;// profit magin

	private boolean parallelComput = false;

	public TestMulticut(int n, int m, int simulationC, boolean parallelComput) {
		this.n = n;
		this.m = m;
		this.simulationCount = simulationC;
		this.parallelComput = parallelComput;
	}

	private void initData() {
		this.capacities = Support.readCapacities(n, m);
		this.demands = Support.readDemands(simulationCount, n);
		this.link_c = Support.readLinkCost(n, m);
		this.profit_m = Support.readProfitMargin(n, m);
	}

	public void detectCriticalS() {
		this.initData();

		CriticalScenarioDetect csd = new CriticalScenarioDetect(n, m,
				simulationCount, demands, capacities);
		List<Double> ratios = csd.maxFlow();
		int above8 = 0;
		int above9 = 0;
		for (double r : ratios) {
			System.out.println(r);

			if (r > 0.85)
				above8++;
			if (r > 0.99)
				above9++;
		}
		System.out.println(above8);
		System.out.println(above9);
	}

	public void Lshape() {
		this.initData();

		double LB = Double.NEGATIVE_INFINITY;
		double UB = Double.MAX_VALUE;
		double gap = Double.MAX_VALUE;

		BigDecimal pro = BigDecimal.ONE.divide(
				BigDecimal.valueOf(simulationCount), 5,
				BigDecimal.ROUND_HALF_EVEN);

		try {
			IloCplex master = new IloCplex();
			master.setOut(new NullOutputStream());
			IloNumVar[][] Y = new IloNumVar[n][m];
			int[][] YValue = new int[n][m];
			IloNumVar[] thetas = master.numVarArray(simulationCount, -1000000,
					1000000);
			double[] thetasValue = new double[simulationCount];
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < m; j++) {
					Y[i][j] = master.boolVar();
					YValue[i][j] = 1;// init Y value
				}
			}
			for (int g = 0; g < simulationCount; g++)
				thetasValue[g] = -100000;// init theta value

			IloLinearNumExpr masterObj = master.linearNumExpr();
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < m; j++) {
					masterObj.addTerm(link_c[i][j], Y[i][j]);
				}
			}
			for (int g = 0; g < simulationCount; g++) {
				masterObj.addTerm(pro.doubleValue(), thetas[g]);
			}
			master.addMinimize(masterObj);

			long start = System.currentTimeMillis();

			int[][] currentOptSolution = Lang.clone(YValue);
			double currentOptVal = 0.0;
			// int delta = n / 4;

			Map<IloRange, Double> dualValues = new HashMap<>();
			List<IloRange> ranges = new ArrayList<IloRange>();

			// // @enhance warm start
			// WarmStartResult wsr = new WarmStartResult();
			// wsr.runWarmStart();
			// currentOptVal = wsr.upperBound;
			// currentOptSolution = Lang.clone(wsr.initSolution);
			// // pending

			int iter = 1;
			while (true) {
				boolean eliminated = true;
				if (iter > 1) {// solve master problem

					// @enhance dynamically adjust tolerance
					if (gap > 0.05)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.08);
					else if (gap > 0.03)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.05);
					else if (gap > 0.02)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.03);
					else if (gap > 0.01)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.02);
					else if (gap > 0.008)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.008);
					else if (gap > 0.007)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.007);
					else if (gap > 0.005)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.005);
					else if (gap > 0.003)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.003);
					else if (gap > 0.002)
						master.setParam(IloCplex.DoubleParam.EpGap, 0.002);
					else
						master.setParam(IloCplex.DoubleParam.EpGap, 0.001);
					// end enhance

					// // @enhance provide initial solution for each iteration
					// for (int i = 0; i < n; i++) {
					// double[] init = new double[m];
					// for (int j = 0; j < m; j++)
					// init[j] = currentOptSolution[i][j];
					// master.addMIPStart(Y[i], init);
					// }
					// // end enhance

					// @enhance warm start with small scale equivalent
					// deterministic problem
					long edpStart = System.currentTimeMillis();
					int smallSize = 50;
					DesignMIP dmip = new DesignMIP(n, m, smallSize, 500);
					int[][] initY = dmip.MIPSolve();
					for (int i = 0; i < n; i++) {
						double[] col = new double[m];
						for (int j = 0; j < m; j++)
							col[j] = initY[i][j];
						master.addMIPStart(Y[i], col);
					}
					long edpEnd = System.currentTimeMillis();
					System.out.println("small scal EDP time = "
							+ (edpEnd - edpStart));

					// @enhance warm start with relaxation
					// long relaxStart = System.currentTimeMillis();
					// double[][] relaxedInit = new double[n][m];
					// IloConversion[][] conversions = new IloConversion[n][m];
					// for (int i = 0; i < n; i++) {
					// for (int j = 0; j < m; j++) {
					// conversions[i][j] = master.conversion(Y[i][j],
					// IloNumVarType.Float);
					// master.add(conversions[i][j]);
					// }
					// }
					// if (master.solve()) {
					// for (int i = 0; i < n; i++)
					// relaxedInit[i] = master.getValues(Y[i]);
					// }
					// for (int i = 0; i < n; i++) {
					// for (int j = 0; j < m; j++) {
					// master.remove(conversions[i][j]);
					// }
					// }
					//
					// for (int i = 0; i < n; i++) {
					// master.addMIPStart(Y[i], relaxedInit[i]);
					// }
					// long relaxend = System.currentTimeMillis();
					// System.out.println("relax time = "
					// + (relaxend - relaxStart));

					// // @enhance: trust region
					// IloRange trustRegion = null;
					// IloLinearNumExpr trustExpr = master.linearNumExpr();
					// int onecount = 0;
					// for (int i = 0; i < n; i++) {
					// for (int j = 0; j < m; j++) {
					// trustExpr.addTerm(1, Y[i][j]);
					// if (currentOptSolution[i][j] == 1)
					// onecount++;
					// }
					// }
					// trustRegion = master.addLe(trustExpr, delta + onecount);
					// // end enhance

					long masterStart = System.currentTimeMillis();

					if (master.solve()) {
						long masterend = System.currentTimeMillis();
						// System.out.println("master time = "
						// + (masterend - masterStart));
						double currentMasterObjV = master.getObjValue();
						if (currentMasterObjV > LB)
							LB = currentMasterObjV;
						else
							eliminated = false;

						for (int i = 0; i < n; i++)
							for (int j = 0; j < m; j++)
								YValue[i][j] = (int) master.getValue(Y[i][j]);
						thetasValue = master.getValues(thetas);
					}

					// if (trustRegion != null)
					// master.delete(trustRegion);
				}

				double estimatedValue = 0.0;
				double subObjValue = 0.0;
				for (int i = 0; i < n; i++) {
					for (int j = 0; j < m; j++) {
						if (YValue[i][j] == 1)
							estimatedValue += link_c[i][j];
					}
				}

				// handle each scenario
				long subStart = System.currentTimeMillis();
				boolean hasMoreCut = false;

				if (parallelComput) {
					// @enhance parallel computing
					Executer executer = new Executer(5, true);
					for (int k = 0; k < simulationCount; k++) {
						executer.fork(new DualSolver(n, m, demands[k], YValue,
								profit_m, capacities, k));
					}
					List<DualResult> results = executer.join();
					executer.shutdown();
					for (DualResult dr : results) {
						if (dr != null) {
							subObjValue += dr.getObjValue();
							if (dr.getObjValue() > thetasValue[dr
									.getTaskIndex()]) {
								hasMoreCut = true;
								BigDecimal[][] Et = dr.getEt();
								BigDecimal et = dr.getEtm();
								IloLinearNumExpr cut = master.linearNumExpr();
								for (int i = 0; i < n; i++) {
									for (int j = 0; j < m; j++) {
										double coeff = Et[i][j].doubleValue();
										cut.addTerm(coeff, Y[i][j]);
									}
								}
								cut.addTerm(1.0, thetas[dr.getTaskIndex()]);
								// add optimality cut
								double rhs = et.doubleValue();
								IloRange range = master.addGe(cut, rhs);
								ranges.add(range);
							}
						}
					}
					// enhance parallel computing end
				} else {
					/* single thread computing begin */
					for (int k = 0; k < simulationCount; k++) {
						IloCplex dual = new IloCplex();
						dual.setOut(new NullOutputStream());
						IloNumVar[] phi = dual.numVarArray(m,
								Double.NEGATIVE_INFINITY, 0.0);
						IloNumVar[] sigma = dual.numVarArray(n,
								Double.NEGATIVE_INFINITY, 0.0);
						IloNumVar[][] gamma = new IloNumVar[n][m];
						for (int i = 0; i < n; i++)
							for (int j = 0; j < m; j++)
								gamma[i][j] = dual.numVar(
										Double.NEGATIVE_INFINITY, 0.0);

						IloLinearNumExpr dualObj = dual.linearNumExpr();
						for (int j = 0; j < m; j++)
							dualObj.addTerm(capacities[j], phi[j]);
						for (int i = 0; i < n; i++)
							dualObj.addTerm(demands[k][i], sigma[i]);
						for (int i = 0; i < n; i++) {
							for (int j = 0; j < m; j++) {
								if (YValue[i][j] == 1)
									dualObj.addTerm(capacities[j], gamma[i][j]);

								// constraint
								dual.addLe(
										dual.sum(phi[j], sigma[i], gamma[i][j]),
										-profit_m[i][j]);
							}
						}
						dual.addMaximize(dualObj);

						if (dual.solve()) {
							double objv = dual.getObjValue();
							subObjValue += objv;
							double[] sigmaVal = dual.getValues(sigma);
							double[] phiVal = dual.getValues(phi);
							double[][] gammaVal = new double[n][m];
							for (int i = 0; i < n; i++)
								gammaVal[i] = dual.getValues(gamma[i]);

							if (objv > thetasValue[k]) {
								hasMoreCut = true;
								BigDecimal et = BigDecimal.ZERO;
								BigDecimal[][] Et = new BigDecimal[n][m];
								for (int i = 0; i < n; i++) {
									for (int j = 0; j < m; j++)
										Et[i][j] = BigDecimal.ZERO;
								}

								for (int i = 0; i < n; i++) {
									et = et.add(BigDecimal.valueOf(sigmaVal[i]
											* demands[k][i]));
								}
								for (int j = 0; j < m; j++) {
									et = et.add(BigDecimal.valueOf(phiVal[j]
											* capacities[j]));
								}
								for (int i = 0; i < n; i++) {
									for (int j = 0; j < m; j++) {
										Et[i][j] = Et[i][j].add(BigDecimal
												.valueOf(gammaVal[i][j]
														* (-capacities[j])));
									}
								}

								IloLinearNumExpr cut = master.linearNumExpr();
								for (int i = 0; i < n; i++) {
									for (int j = 0; j < m; j++) {
										double coeff = Et[i][j].doubleValue();
										cut.addTerm(coeff, Y[i][j]);
									}
								}
								cut.addTerm(1.0, thetas[k]);
								// add optimality cut
								double rhs = et.doubleValue();
								IloRange range = master.addGe(cut, rhs);
								ranges.add(range);
							}
						}
						dual.end();
					}// end for all scenarios
					/* single thread computing end */
				}

				long subEnd = System.currentTimeMillis();
				// System.out.println("subproblem time=" + (subEnd - subStart));

				estimatedValue += BigDecimal.valueOf(subObjValue).multiply(pro)
						.doubleValue();

				// @enhance
				if (estimatedValue < currentOptVal) {
					currentOptVal = estimatedValue;
					currentOptSolution = Lang.clone(YValue);
				}
				// else {
				// // @enhance eliminate current solution
				// if (iter > 3) {
				// IloLinearNumExpr cutExp = master.linearNumExpr();
				// int countOne = 0;
				// for (int i = 0; i < n; i++) {
				// for (int j = 0; j < m; j++) {
				// if (YValue[i][j] == 1) {
				// countOne++;
				// cutExp.addTerm(-1, Y[i][j]);
				// } else
				// cutExp.addTerm(1, Y[i][j]);
				// }
				// }
				// master.addGe(cutExp, 1 - countOne);
				// }
				// }

				System.out.println("estimatedValue=" + estimatedValue);
				if (estimatedValue < UB)
					UB = estimatedValue;
				else
					eliminated = false;

				// if (eliminated) {
				// // @enhance eliminate current solution
				// IloLinearNumExpr cutExp = master.linearNumExpr();
				// int countOne = 0;
				// for (int i = 0; i < n; i++) {
				// for (int j = 0; j < m; j++) {
				// if (YValue[i][j] == 1) {
				// countOne++;
				// cutExp.addTerm(-1, Y[i][j]);
				// } else
				// cutExp.addTerm(1, Y[i][j]);
				// }
				// }
				// master.addGe(cutExp, 1 - countOne);
				// }

				if (iter > 1)
					gap = BigDecimal
							.valueOf(UB - LB)
							.divide(BigDecimal.valueOf(-LB), 5,
									BigDecimal.ROUND_HALF_EVEN).doubleValue();

				System.out.println("iter=" + iter + "\t UB=" + UB + "\t LB="
						+ LB + "\t gap=" + gap);
				long end = System.currentTimeMillis();
				long time = end - start;

				System.out.println("total time = " + time);

				if (gap < 0.001 || time > 3600000 || !hasMoreCut) {
					StringBuilder builder = new StringBuilder("multicut; DAG");
					if (parallelComput)
						builder.append("parallel comput; ");
					builder.append("n=" + n + ",m=" + m + ",M="
							+ simulationCount + "\r\n");

					for (int i = 0; i < n; i++) {
						for (int j = 0; j < m; j++) {
							if (YValue[i][j] == 1) {
								System.out.print("(" + i + "," + j + ")\t");
								builder.append("(" + i + "," + j + ")\t");
							}
						}
						System.out.println();
						builder.append("\r\n");
					}
					System.out.println("Time consumed = " + time);
					double value = BigDecimal
							.valueOf(LB + UB)
							.divide(BigDecimal.valueOf(2), 5,
									BigDecimal.ROUND_HALF_EVEN).doubleValue();
					builder.append("Time consumed = " + time + "\t gap=" + gap
							+ "\t value=" + value);
					builder.append("\r\n");

					try {
						Writer writer = new BufferedWriter(
								new FileWriter(
										Lang.projectRoot()
												+ "\\data\\flexibilityDesign\\result\\result.txt",
										true));
						writer.write(builder.toString());
						writer.close();
					} catch (IOException e) {
						// TODO Auto-generated catch block
						e.printStackTrace();
					}

					break;
				} else
					iter++;
			}
			master.end();

		} catch (Exception e) {
			e.printStackTrace();
		}
	}

	private class WarmStartResult {
		private int warmCount = 5;// warm run times
		private int warmSize = 10;// instances for each warm run
		private int[][] randomIndexs = new int[warmCount][warmSize];

		double upperBound = Double.MAX_VALUE;
		int[][] initSolution = new int[n][m];

		private void generateRandomIndex() {
			int max = demands.length;
			Random random = new Random();
			for (int i = 0; i < warmCount; i++) {
				for (int j = 0; j < warmSize; j++) {
					randomIndexs[i][j] = random.nextInt(max);
				}
			}
		}

		void runWarmStart() {
			System.out.println("@ begin warm start");
			generateRandomIndex();

			double pro = BigDecimal.ONE.divide(BigDecimal.valueOf(warmSize), 5,
					BigDecimal.ROUND_HALF_EVEN).doubleValue();

			for (int wi = 0; wi < warmCount; wi++) {
				try {
					IloCplex master = new IloCplex();
					master.setOut(new NullOutputStream());
					IloNumVar[][] Y = new IloNumVar[n][m];
					for (int i = 0; i < n; i++)
						for (int j = 0; j < m; j++)
							Y[i][j] = master.boolVar();

					IloNumVar[][][] x = new IloNumVar[warmSize][n][m];
					for (int sm = 0; sm < warmSize; sm++) {
						for (int i = 0; i < n; i++) {
							for (int j = 0; j < m; j++) {
								x[sm][i][j] = master.numVar(0.0,
										Double.MAX_VALUE);
							}
						}
					}
					IloLinearNumExpr obj = master.linearNumExpr();
					for (int i = 0; i < n; i++)
						for (int j = 0; j < m; j++)
							obj.addTerm(link_c[i][j], Y[i][j]);
					for (int sm = 0; sm < warmSize; sm++)
						for (int i = 0; i < n; i++)
							for (int j = 0; j < m; j++)
								obj.addTerm(-profit_m[i][j] * pro, x[sm][i][j]);
					master.addMinimize(obj);

					for (int sm = 0; sm < warmSize; sm++) {
						for (int j = 0; j < m; j++) {
							IloLinearNumExpr ca = master.linearNumExpr();
							for (int i = 0; i < n; i++) {
								ca.addTerm(1.0, x[sm][i][j]);
							}
							master.addLe(ca, capacities[j]);
						}
						for (int i = 0; i < n; i++) {
							IloLinearNumExpr de = master.linearNumExpr();
							for (int j = 0; j < m; j++) {
								de.addTerm(1.0, x[sm][i][j]);
							}
							master.addLe(de, demands[randomIndexs[wi][sm]][i]);
						}
						for (int i = 0; i < n; i++)
							for (int j = 0; j < m; j++) {
								IloLinearNumExpr cab = master.linearNumExpr();
								cab.addTerm(Y[i][j], capacities[j]);
								master.addLe(x[sm][i][j], cab);
							}
					}

					double[][] init = new double[n][m];
					for (int i = 0; i < n; i++)
						for (int j = 0; j < m; j++)
							init[i][j] = 1.0;
					for (int i = 0; i < n; i++)
						master.addMIPStart(Y[i], init[i]);

					if (master.solve()) {
						double objV = master.getObjValue();
						if (objV < upperBound) {
							upperBound = objV;
							for (int i = 0; i < n; i++) {
								for (int j = 0; j < m; j++) {
									initSolution[i][j] = (int) master
											.getValue(Y[i][j]);
								}
							}
						}

					}
					master.end();
				} catch (Exception e) {
					e.printStackTrace();
				}
			}
			System.out.println("@ end warm start");
		}
	}
}
